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ABSTRACT 

The interaction of accretion disks with the magnetospheres of young 
stars can produce X-winds and funnel flows. With the assumption of 
axial symmetry and steady state flow, the problem can be formulated 
in terms of quantities that are conserved along streamlines, such as the 
Bernoulli integral (BI), plus a partial differential equation (PDE), called 
the Grad-Shafranov equation (GSE), that governs the distribution of 
streamlines in the meridional plane. The GSE plus BI yields a PDE of 
mixed type, elliptic before critical surfaces where the flow speed equals 
certain characteristic wave speeds are crossed and hyperbolic afterward. 
The computational difficulties are exacerbated by the locations of the 
critical surfaces not being known in advance. To overcome these obsta- 
cles, we consider a variational principle by which the GSE can be attacked 
by extremizing an action integral, with all other conserved quantities of 
the problem explicitly included as part of the overall formulation. To 
simplify actual applications we adopt the cold limit of a negligibly small 
ratio of the sound speed to the speed of Keplerian rotation in the disk 
where the X-wind is launched. We also ignore the obstructing effects 
of any magnetic fields that might thread a disk approximated to be in- 
finitesimally thin. We then introduce trial functions with adjustable 
coefficients to minimize the variations that give the GSE. We tabulate 
the resulting coefficients so that other workers can have analytic forms 
to reconstruct X-wind solutions for various astronomical, cosmochemical, 
and meteoritical applications. 



Subject headings: stars: pre-main-sequence; winds; ISM: accretion disks; 
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1. Introduction 

Accretion, disks, and jets are ubiquitous in astrophysics (see, e.g., Blandford 
& Rees 1992). A consensus has been reached that an extra needed ingredient to 
obtain outflow from inflow is the presence of strong magnetic fields that thread 
a disk conventionally assumed to be rotating at Keplerian speeds about a central 
gravitating object, taken in this paper to be a newly born star. Differences come 
in ascribing the origin of the magnetic fields in the disk itself or in the central star 
(K6nigl & Pudritz 2000, Shu et al. 2000). 

Disk winds have been extensively studied, both analytically via the assumptions 
of self-similarity in 2-D space for axisymmetric, time-independent flows (e.g., Bland- 
ford & Payne 1982; Contopoulos & Lovelace 1994) or by taking advantage of arbi- 
trary variations of the gas pressure (e.g., Tsinganos & Trussoni 1991) or by studying 
the asymptotic properties of the collimation (Heyvaerts & Norman 1997); and nu- 
merically by finite-element methods attacking the axisymmetric, time- independent, 
Grad-Shafranov equation (e.g., in the relativistic regime by Camenzind 1987) or by 
finite-difference treatments of the time-dependent equations of ideal MHD in 2- and 
3-D (e.g., Uchida & Shibata 1986, Pudritz et al. 2006). For a review of these types 
of calculations, see Ferreira (2004). 

The most highly developed semi-analytic theory for the second viewpoint is 
called X-wind theory, in which fast jets arising in young stellar objects (YSOs) owe 
their existence to the interaction of the accretion disk with the magnetosphere of 
the central star. The interaction of accretion disks with strongly magnetized central 
stars has also been studied numerically (e.g., Goodson, Bohm, & Winglee 1999; Long 
et al. 2005, Ustyugova et al. 2006). Although both funnel flows and X-like winds 
have been found, they have yet to appear simultaneously in numerical simulations, 
probably because the numerical calculations have not yet proceeded to steady state 
where the condition of disk-locking applies (Shu et al. 1994). Pure X-wind theory 
assumes for simplicity that the disk itself is unmagnetized, in fact, all that is needed 
for the theory to work is for open field lines to be concentrated in a narrow annulus 
near the inner edge of an accretion disk. 

Recently, Bacciotti et al. (2002) and Coffey et al. (2004) identified jet rotation 
in four T Tau systems, DG Tau, TH 28, RW Aur, and LkHa 321, of an amount too 
large to be compatible with X-winds, but consistent with launching from disks at 
radii of 0.5-2 AU. Later, Cabrit et al. (2006) showed from mm-wave radio measure- 
ments that the disk rotation in RW Aur is actually in the opposite sense to that 
deduced for the jet from optical lines. Moreover, Pety et al. (2006) find that HH 
30, which is observed nearly edge-on and therefore should have had the clearest sig- 
nature for jet rotation, showed no evidence for outflow rotation at mm wavelengths, 
a conclusion reinforced by optical and ultraviolet observations of the HH 30 jet by 
Coffey et al. (2007). While the positive results remain for the three other systems, 



- 3 - 



the case of HH 30, where longitudinal velocities occur in the direction transverse to 
the line of sight, suggests that the slight line asymmetries in the other cases may 
be more associated with unequal jumps in the velocity of shocked, high-speed jets, 
than to the rotation of collimated outflows. 

In contrast, no one has proposed any explanation other than X-winds for the 
correlated inflow-outflow signatures seen in SU Aur by Giampapa et al. (1993) 
and Johns & Basri (1995). Apart from SZ 68 (Johns-Krull & Hatzes 1997), we are 
unaware of any other T Tau star that shows a tilted-dipole magnetic-field geometry, 
and it could be that the dipole component is small on the surface of most T Tau 
stars (Johns-Krull 2007). Fortunately, Mohanty & Shu (2007) show that while funnel 
flows are sensitive to the detailed assumptions made concerning multi-pole structure 
on the surfaces of the central stars, the properties of the X-wind depend mostly only 
on the amount of trapped flux in the X-region (see also the observational evidence 
relating to this point collected by Johns-Krull & Gafford 2002). 

Recent calculations show that YSOs are unlikely to lose enough magnetic flux in 
the process of gravitational collapse to make the level of magnetization ignorable in 
the resultant circumstellar accretion disks (Galli et al. 2006; Shu et al. 2006, 2007). 
Indeed, the disks are sufficiently magnetized in many cases that, in quasi-steady 
state, they rotate at sub-Keplerian rates until the the inner disk-edge is reached. 
Thus, there are open questions of how much of the trapped flux near the inner edge 
is to be attributed to the central star versus the disk, and how such disks reacquire 
near-Keplerian rates of rotation at their inner edges. We ignore these complications 
in the present study of the X-wind phenomenon, but we note that the methods 
introduced here are easily modified to attack the more complex problem when the 
accretion disk interacting with a stellar magnetosphere is itself strongly magnetized. 

The original X-wind model supposed the outflow to occur from the equator of 
a magnetized star spun to breakup by a presence of an accretion disk that abutted 
its surface (Shu et al. 1988). Later, in order to accommodate the slow rotators, 
such as the classical T Tauri stars which are only rotating at one tenth of breakup 
(Vogel & Kuhi 1981, Bouvier et al. 1991, Edwards et al. 1993), Shu et al. (1994a) 
generalized the X-wind picture to include the case of relatively low accretion when 
the magnetosphere of the star would truncate the accretion disk at an inner edge 
before the disk reached the stellar surface (typically a circle of radius 0.2 on the scale 
of Fig.l, where the disk's inner edge is taken to be at w — 1). In a quasi-steady- 
state where most of the mass of the central star is built up by disk accretion, the 
magnetic coupling between the star and the disk regulates the star to corotate at the 
Keplerian frequency at the truncation radius. For a protostar with magnetic dipole 
moment /x*, mass M*, mass accretion rate M D , Ostriker & Shu (1995) estimate this 
radius to be 




(1-1) 



-4 - 



where $dx is an order unity dimensionless number that parameterizes the amount of 
stellar magnetic flux that is trapped in the disk. Inside this radius, matter is chan- 
neled to the star via a funnel flow. The excess angular momentum of the accreting 
material is deposited in the magnetic field in the form of Maxwell torque, and then 
transported back to the disk. The gain of angular momentum and approximate field 
freezing would try to move the footpoint of the funnel-flow field lines outward. 

Exterior to the truncation radius Rx, the equatorial inward drift in the accretion 
disk creates an angle between the stellar magnetic-field lines and the disk normal. If 
approximate field freezing holds as the accretion proceeds, a fraction of the field lines 
will develop an angle larger than 30°, when matter frozen to this flux tube becomes 
unstable to magnetocentrifugal fling (Blandford & Payne 1982). These field lines 
are thus responsible for driving a magnetohydrodynamic (MHD) wind from the 
disk. Since the wind removes angular momentum from the disk, the footpoints of 
those field lines in the disk will try to migrate inward. The radially inward press 
of the footpoints of the wind field lines footpoints and the radially outward press 
of the footpoints of the funnel-flow field lines create a magnetic X-configuration 
that distinguishes the model from similar variants in the literature (the historical 
choice of the name from the X-point of the equivalent gravitational potential in the 
co-rotating frame is common to many models). In quasi-steady state where radial 
advection into the X is balanced by the resistive diffusion of field lines out of the 
X, Shu et al. (1994b) estimated that enough stellar flux could be trapped in a 
small X-region near the inner disk edge to have large dynamical effects, namely, the 
truncation of the disk by a funnel flow out of the disk plane accompanied by an 
X-wind that carries away most of the excess angular momentum transported into 
the X-region. 

Apart from the original numerical estimates, there are reasons to suppose that 
if turbulent resistivity is the source of the diffusion across magnetic field lines (Shu et 
al. 2007), then the fractional size of the X-region in units of Rx is given by the ratio 
of sound speed at the surface of the disk where the X-wind is launched to the local 
Keplerian speed at R X - For the inner disk of a classical T Tauri star, the thermal 
sound speed is a ~ 5 km/s while the Keplerian speed at Rx is vk ~ 100 km/s 
(Najita et al. 2007); thus, the ratio e is a small number ~ 0.05. In an asymptotic 
analysis where e is taken to — > 0, the X-wind tied to the trapped field lines in the 
X-region would emerge from virtually a single point in the meridional plane with a 
fan-like geometry. Seen by an observer rotating at the Keplertian angular frequency 
of Rx, gas flows along streamlines that coincide with field lines if field freezing is 
assumed, and both patterns of streamlines and field lines remain stationary in the 
corotating frame. 

Viewed in this fashion, the overall problem can be broken into smaller pieces and 
tackled separately. Using a formulation with a precedent in the work of Lovelace 
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et al. (1986), Shu et al. (1988, 1994b) wrote down the mathematical equations 
that describe a steady state axisymmetric flow in the corotating frame (Grad & 
Rubin 1958, Shafranov 1966). By the method of matched asymptotic expansions, 
Shu et al. demonstrated the existence of an inner solution in the X-region where 
the flow makes a sonic transition, and an outer solution where the sound speed is 
formally taken to zero and the X-region shrinks down to a point. Najita & Shu 
(1994) computed numerically the portion of the X-wind in which the fluid velocity 
is sub-Alfvenic, and the governing equation is elliptic. Ostriker & Shu (1995) solved 
the problem of the funnel flow and the field configuration in the dead zone (in which 
the field lines do not depart sufficiently from disk normal to load any matter) in an 
approximation that treated the accretion flow onto the star as highly sub-Alfvenic. 
Shu et al. (1995) constructed asymptotic solutions that describe the logarithmically 
slow, far-wind collimation into jets. The free boundaries between various parts of the 
problem (funnel flow, dead zone, and X-wind) were determined by pressure balance 
on either side. 

In this paper, we wish to address the X-wind part of the overall problem. In 
order for the X-wind to accelerate from rest to supersonic speeds, it must smoothly 
pass three surfaces on which the flow velocity is equal to the slow MHD, Alfven, and 
fast MHD velocity, respectively (Heinemann & Olbert 1978, Sakurai 1985). These 
critical surfaces manifest themselves as singularities in the governing equation (see 
Weber & Davis 1967), and thus need to be handled analytically. In an axisymmet- 
ric problem, if the shapes of the streamlines are known in the meridional plane, 
the conserved quantities of the problem (mass to flux loading, angular momentum 
flux including that carried in the Maxwell stress, and Bernoulli's integral along a 
streamline) suffice to give a completely analytic solution, including the locations 
and conditions required to cross the critical surfaces smoothly. Unfortunately, the 
streamline distribution in the meridional plane is not known a priori but must be 
obtained, in principle, from a solution of the Grad-Shafranov equation (GSE). The 
spatial location of the critical surfaces are part of the overall solution of the GSE; 
indeed, they characterize the regions where this PDE is elliptic or hyperbolic. The 
mixed character of the GSE makes a direct numerical attack extremely difficult 
when self-similarity does not apply, perhaps the hardest problem in the mathemat- 
ical theory of nonlinear PDEs of second-order (see Garabedian 1986). The current 
work sidesteps the mathematical solution of the GSE as a nonlinear PDE of second 
order, and approaches it instead as a much more amenable problem in variational 
calculus. 

The rest of this paper is organized as following. In §2 we review the basic for- 
mulation in terms of a stream function and an Alfven discriminant that yields the 
partial-differential equations - the so-called Grad-Shafranov and Bernoulli equations 
- that govern a steady, axisymmetric, X-wind flow. In §3 we write down an action 
whose variations with respect the stream function ip{w, z) and the Alfven discrim- 



- 6 - 



inant A(to, z) yield, respectively, the Grad-Shafranov equation and the Bernoulli 
equation when the action reaches an extremum. We also perform a transformation 
where we replace the vertical coordinate z in cylindrical coordinates (w, </?, z) by ip. 
In §4, we show how to incorporate boundary conditions into the problem, as well as 
how to take advantage of the fact that analytic forms are known for the solutions 
in the near-neighborhood of the X-point and in the asymptotic regime far from the 
X-point (Shu et al. 1994b, 1995). In §5, we outline a practical implementation 
of the principle of extremal action, making use of only variations of - or, more 
precisely, of z(w, ip) in our actual working space - as the substitute to attacking the 
Grad-Shafranov equation, while we solve Bernoulli's equation directly for reasons 
that are expounded upon in this section. In §6, we present numerical results for 
three specific cases of mass-loading onto wind flux-tubes, finding good agreement 
with previous approximate solutions obtained by Shang (1998) that have been used 
for many different astrophysical applications (e.g., Shang et al. 1998, 2002, 2004). 
In §7, we summarize the recipes needed to convert the numerical solutions of §6 
into practical dimensional models. We then offer our conclusions and suggestions 
for needed future research. 



2. Basic Equations 

From the fundamental parameters of the problem, we may construct units of 
length, time, and density as Rx, an d M w /4nR x Qxi respectively. By assuming 
axisymmetry and stationarity in a frame that is rotating with angular velocity fix, 
we may write down the dimensionless governing equations in the above units. 

V • (pu) = 0, (2-la) 

V(^|u| 2 ) +(2e z + Vxu) xu = --Vp-W cff + i(VxB) x B,(2-lb) 

V 2 / P P 

B x u = 0, (2-lc) 

V • B = 0, (2-ld) 

where e = a/Rx^lx is the sound speed measured in units of Keplerian velocity at 
the X-point, and is assumed to be a small parameter of the problem. The effective 
potential in the corotating frame is defined as 

V cS = 1 + - (2-2) 

Vw 2 + z* 2 2 K J 



Here we have added a constant term to the effective potential so that its numerical 
value is zero at the X-point. 
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2.1. Constants of Motion 

The continuity equation (12- lap is satisfied identically if we define the poloidal 
velocity through a stream function (Shu et al. 1988, 1994a): 

1 dip 1 dip 
pu™ = — -5-, pu z = ^— • (2-3) 

For steady state axisymmetric flow in the corotation frame, the field freezing con- 
dition (12- left demands that the magnetic field and the velocity are related by (see, 
e.g., Mestel 1968) 

B = (3pu. (2-4) 

With this identification, the continuity equation (12- lap and the absence of magnetic 
monopoles (12-ldl) imply u • V/3 = 0. In terms of the stream function, this means (3 
is constant along each streamline, or (3 = /3(ip). 

The Euler equation describes momentum and energy balance in three spatial 
dimensions. If we take the component along the fluid velocity by taking the inner 
product of (12-lbl) with u, we obtain the Bernoulli's equation (BE) along streamlines 



u- VH = where H = ^|u| 2 + e 2 hip + V cfr . (2-5) 

In other words, H = H(ip), and the energy per unit mass of an isothermal gas, 
including its specific enthalpy, is conserved along a streamline in the corotating frame 
where the flow occurs parallel to B. Similarly, if we take the toroidal component 
of the Euler equation ( 12-lbl) . we obtain a third conserved quantity along stream 
lines, the angular momentum of the gas allowing for that part carried away by the 
Maxwell torque of the field: 

J = zu 2 + w(l - P 2 p)u !fi = J(ip). (2-6) 

As we shall see, the determination of the conserved quantities, H(ip)and J(ip), is 
achieved by demanding that the X-wind crosses the slow MHD and fast MHD sur- 
faces smoothly. The loading of mass onto flux, which is governed by /3(ip) is freely 
specifiable within certain limits to be detailed below. 

The last component of the Euler equation describes momentum balance in the 
direction perpendicular to the poloidal field lines. It is the famous Grad-Shafranov 
equation (Heinemann & Olbert 1978, Sakurai 1985): 

, AT7 ,s If J 1 \J , P 2 '{V eS + e 2 In [6 2 h/((3 2 - w 2 A)]} e 2 h'/h 

(2-7) 

where we have rescaled Bernoulli's function as 



H = -e 2 ln(e 2 /i), 



(2-8) 
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so that h remains an order unity quantity in our calculation. Here A is the Alfven 
discriminant defined by 

A=*£f±, (2-9) 

where 

Ml s £ = JL (2 . 10) 

is the Alfven Mach number. Hence A is positive when the total velocity is less than 
the Alfven speed, and negative when the total velocity is larger than the Alfven 
speed. From the form of the GSE (12-71) . we see that the conserved angular momen- 
tum flux J is not freely specifiable. It is determined by the condition of smooth 
Alfven transition. In order for the solution to remain continuous and differentiable, 
one must impose 

J = w 2 whenever A = 0. (2-11) 

The elimination of p in the equations in favor of A is based on numerical considera- 
tions, since p will in general vary by many orders of magnitude, while A only varies 
moderately. As we argued in the previous section, the sound speed e is likely to be 
small. In terms of these variables, the BE takes the form 

|2 , 1 ( J X 2w 2 {V cS + e 2 \n[e 2 h/((3 2 -zu 2 A)]} 



™ + tA^-- 1 ) + V-^T' — = °- (2 ' 12) 



2.2. The Cold Limit 

With A implicitly defined in the BE ( 15121) . the GSE is a PDE of mixed type, 
which demands different numerical methods in different regions (see Heinemann & 
Olbert 1978 and Appendix [A]) . There are three relevant signal speeds (which we 
term sonic, slow, and fast in the appendix) involved in an MHD flow (see Jackson 
1975 or Shu 1992). The loci where the poloidal fluid speed equals those signal 
speeds separate the flow into four regions. As the poloidal velocity exceeds the 
sonic speed, the governing GSE changes from elliptic to hyperbolic. A wise strategy 
might start with the search of appropriate boundary conditions in the disk where 
u 2 = 0, and at the sonic surface (whose location is still undetermined), followed 
by a standard scheme (e.g., relaxation) to obtain the interior solution. Beyond the 
sonic surface, the GSE becomes hyperbolic. The boundary condition on the sonic 
surface now serves as the initial condition, which we use to integrate forward along 
characteristics toward the slow surface. We then follow similar procedures to obtain 
solutions from the slow surface to the fast surface, and beyond. 

A significant simplification can be achieved when the sound speed is negligible, 
as in the outer problem of the X-wind (see §4 of Shu et al. 1994b). The governing 
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equations are treated as power series expansion in e. The leading term in the GSE 



Notice the lowest order term in H vanishes independent of the form of h. In this 
limit, both the sonic speed and the slow speed reduce to zero, and the first elliptic 
and hyperbolic parts of the flow shrink down to the X-point. We are thus spared 
the vicissitudes of this portion of the problem. Once the fluid leaves the X-region 
(with poloidal velocity greater than the slow speed), it proceeds to the fast surface, 
where the governing equation becomes hyperbolic. 

Najita & Shu (1994) solved the GSE in the sub-Alfvenic region. By introducing 
a generalized coordinate system, they were able to map the location of the Alfven 
surface to a known location, and determine the functional form of ft (if}) based on 
the position and shape of the Alfven surface. Their numerical scheme to find ft (if)) 
by iteration encountered a systematic "drift problem", however, and an artificial 
"Alfven seam" was invented to cope with this difficulty. 

In a later treatment by Shang (1998), ft (if)) was specified in advance, limited in 
its functional form by considerations of how the gas exits the X-region, an analysis 
that we repeat in §4.2 (see also §5). The GSE was not solved as a PDE, but rather 
as an error estimator in a Weber-Davis type of analysis, where if) as a trial function 
of spatial location is obtained by interpolating between the known analytic forms 
in the X-point neighborhood (see §4.1) and at asymptotic infinity (see §4.3). The 
interpolation formula has a number of degrees of freedom, which are adjusted to 
give "least error" in some sense when the trial solution for if> is substituted back 
into the GSE. The rest of the problem, including the constraints of the conserved 
quantities and smooth passage through the Alfven and fast surfaces, are performed 
exactly. She verified the result derived by Goldreich and Julian (1970) that passage 
through the Alfven surface is automatic in such a scheme if one has guaranteed it 
through the fast surface. In fact, §5.2 demonstrates the falsity of the frequent claim 
made otherwise in the literature that J(if>) is set at the Alfven surface; the claim 
holds only if one already has a solution such that the wind passes smoothly through 
the fast surface. 



Based on the above arguments, the X-wind is a fierce mathematical beast, and a 
direct numerical attack is unlikely to subdue it fully. To construct a global solution of 



(EZZD and the BE (EE2) are 




(2-13b) 



(2-13a) 



3. Variational Principle 
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the X-wind that accelerates elements of plasma from the disk to super-magnetosonic 
speeds, we must resort to a different approach. Consider the following action written 
down by inspect ion Q 



S 




A\Vij>\' 



-C--iY 

2A\w 2 J 



V eS + e 2 \n[e 2 h/{[3 2 - w 2 A)\ - e 2 
I3 2 - tu 2 A 



>GT£, 



(3-1) 

It is straight forward to demonstrate that variation against ip yields the GSE (12-71) . 
while variation against A gives the BE (12-121) . The challenges of constructing solu- 
tions to a nonlinear PDE of mixed type is now transformed to tuning trial functions 
of ip and A until a local extremum of the action (13- ip is reached. 

To formulate a scheme that is easy to implement numerically, we consider a 
change of independent variables from the usual cylindrical coordinates 

(w,z,tp) -> (w,ip,(f). 

For a given value of ip, the functional form of z(w) determines the shape of the 
given streamline, and A(w) offers information on the velocity distribution along 
that streamline. Written in these new coordinates, and taking the cold limit as 
e — > 0, the action reads 



S 




1+1^ 







W 



2Ad^p\w^ ~ 



+ 



dz V e 



eff 



dip (3 2 - w 2 A 



>zudipdzu 



(3-2) 

Since A only enters the action as a constraint rather than a dynamic variable (i.e., 
its derivative is absent in the action), variation with respect A yields the BE as 
before, but now written in a different set of coordinates, 



(dz_ 



1 + 



fdz_^ 2 



Variation with respect to z gives 

ss z 



2tt 



1 dbz(J_ 



2A di\) \w 2 



1 

^ 2 V 



dz\ 

2 

+ 



'd8z 



dip 
d5z 



1 + 



2w 2 V cS 
((3 2 - w 2 Af 



(3-3) 



dz 
dw 



+ A 



dz 
dip 



dz \ d5z 



dz V e 



efF,z 



dip (3 2 - tu 2 A diP(3 2 - w 2 A 



dw ) dw 



Sz \wdipdw. 



l \n their pioneering development of magnetized stellar winds in a Grad-Shafranov formalism, 
Heinemann & Olbert (1978) noted in passing that the resultant equations could be derived from a 
principle of least action, which differs in detailed form from that used in this paper, but is in the 
same spirit. However, they, and subsequent workers who have made similar observations, did not 
exploit the principle to obtain actual wind solutions. 
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Integrating by parts, we have 



5S? = 2tt 



Ihp) 



-2 



fdz^ 2 
\dw 



-(- 

2 A \w 2 



1 + 



cff 



2 



W 



■0 = 



dip 



VD = \ 



+2vr 



#0 / 
dz wV cS , 



1 + 



9 



d0 /3 2 - zu 2 A dw 



A 



dz \ 



+ 



-i 



.(. 



J 



2A\w 2 



cff 



(3 2 - t?M 



wSzdipdw 



- 



dip J 



dz 
dzu 



w 



5zdipdw. 



Since we specify the boundary condition z = on ip = 0, and z = Z{w) on ip = 1, 
where Z(w) is a known function, we see that <5z vanishes on these two boundaries. 
As we shall see later (see |48|; |49| ) . the solution near the X-point and asymptotically 
can be constructed analytically. Thus 5z also vanishes when w — 1 and w — > oo 
in our variational scheme, and both surface terms vanish in the above expression. 
In order for the action to be stationary against any choice of 5z, the solution must 
satisfy the Euler Lagrange equation, 



5S_ 

Tz 




1 + 



dz 
dw 



dz w dV e s d 
dip (3 2 — w 2 A dz dw 



A 



dz 
dip 



1 




2A 


- 

\w 2 




'dz" 




dw . 



V, 



cff 



(3 2 - w 2 A 



■w 



w 



(3-4) 



Dividing both sides by w, we can simplify the Euler-Lagrange equation f)3-4p to 
obtain 



2lhp\ \dip 



\dw 



--(- 

A\w 2 



1 



J' 



+ 



A 2 \w 2 J 



2w 2 V cS 
(/? 2 - w 2 A) 2 



1 f dz 



+ 



w 2 {(3 2 -w 2 A) 2 w\dip 



-i 



_d_ 

dw 



A 



dz 
dw 



w 



\dip J dip] \dip J 



1 + 



(-) 

\dw 



0. 



(3-5) 



We notice that the coefficient of dAjdip is simply the BE, which vanishes at a 
local extremum of the action. One may easily check that the other terms yield the 
conventional GSE, but written in our new coordinates. 
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4. Boundary Conditions 
4.1. X-point 

With the new coordinates, the computational domain is bounded by if) G [0, 1] 
and V3 G [1, oo). The X-point in these coordinates is a singularity given by w = 1 
for all values of if). Fortunately, we have analytic solutions there. From this point 
onward, we shall work with a scaled Alfven discriminant 



This function has the advantage of remaining finite even when f3 diverges. For a 
given functional form of (3 (which tells us how matter is loaded onto the field lines), 
the Alfven discriminant has the series expansion in w — 1 . 

Xx = 1 + xMi™ - 1) + X 2 (V0O - I) 2 + •••> (4-2) 

Here a subscript X reminds us that this series solution is valid near the X-point. 
In order to match asymptotically onto the outer limit of the inner problem (Shu et 
al. 1994b), the density p = /3~ 2 (1 — m 2 x)~ 1 must diverge as (w — l) -2 near the 
X-point. This requirement translates to Xi = ~2 for all values of if). Similarly, for 
the coordinate z (now a dependent variable), we expand it as, 

z x = Zl (iP)(w - 1) + z 2 (if))(w - l) 2 + .... (4-3) 

Notice that the Jacobian near the X-point is 

yfg = z[(w — 1) as w — > 1. 

which is expected since the entire line of if) G [0, 1] is mapped into a single point, 
and the Jacobian must vanish in this situation. Substituting the series expansion 
into the GSE (13-51) in the transformed coordinates, the lowest order is (w — 1)~ 2 . 



d_ 

dip 



4(i + *2) 



0. 



which has the solution 



1 



zx = tarn?, d = — / (3dip. (4-4) 







If we assume that the upper boundary of the X-wind (if) = 1) near the X-point forms 
an angle of 9x with the x axis, we have 



tan6>x ~ 



1 r 

— Z\(i\) — 1) = tan — / pdif). 
1 K Jo 
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Thus, the integration constant is given by 



K = A (3= f pdt/j. (4-5) 



®x Jo 

Substituting back into the Bernoulli's equation allows us to solve for \2 



V4cos 2 ^-l 
Kp cos z v 

For better numerical accuracy, we carry out the computation to next order. The 
next term in the series expansion of the GSE is 0[(w — l) -1 ]: 

d 2 Q 

— + Q = sintf, 

where <2 = z<i cos 3 Given the boundary condition z = at ip = for all values of 
X, the above equation may be solved to give 

z 2 = i S ec 2 $(gtan$-$). (4-7) 

The integration constant q can be determined by expanding the upper boundary 
near the X-point. Substituting into the BE, we can determine the last term without 
the knowledge of J as 



cos 2 & — 2 + tan^((7tan^ — ■&) yj 4 cos 2 •& — 1 (4 — q sec 2 ■&) , oN 
y-x = — — - H — 4. (4-8 ] 

2K f3 cos 2 cos 2 1 2^/3 cos 2 ■& v ; 



4.2. Specifying Mass Loading and Difficulties with the Boundary Layer 

Since the X-wind is driven magnetocentrifugally, one would naively expect that 
it is bounded away from the polar axis (at least in the immediate vicinity of the 
X-point) by some curve which intersects the disk. In the outer limit of the inner 
problem (see eq. 3.10d of Shu et al. 1994b), the gas pressure p = e 2 p takes the form 

(X~ 2 (4cos 2 $ - 1)~ 1/2 , as a -> oo, 

where tan^ = zjiw — 1). As d x — > tt/3 (which is the critical angle for the last 
matter carrying streamline), the pressure diverges unless 

P oc (4 cos 2 tf- 1)~ 1/2 , as V -> 1. 

Substitute this functional form into the lowest order equation ( 14-41) . finite magnetic 
field and pressure on the last streamline demands 



g 



f3oc (1- V)~ 1/3 , as V -> 1. 



(4-9) 
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The divergence of (5 should not come as a surprise. Recall that the last stream- 
line is defined to be the boundary between the X-wind and the dead zone. To ensure 
analyticity across this boundary, we must have p — > as tp — > 1. Now since neither 
the magnetic field nor the velocity become singular, we must take (3 — > oo on that 
last streamline, so that the product f3 2 p 2 = B 2 /u 2 remains finite. With this limit in 
place, we see that the rescaled Alfven discriminant 

1 ~ 1/P 2 P 1 

X 9 ^ 9 

W ZD 

remains positive for all points along the last streamline. In other words, the flow on 
the last streamline is always sub-Alfvenic, since the Alfven speed is infinite there. 
This behavior of the last streamline requires a double limiting procedure if we were 
to accurately construct the asymptotic solution on that interface. We speculate that 
this difficulty is an indication that the last streamline needs to be treated as a bound- 
ary layer. This speculation is reinforced by the fact that the last X-wind streamline 
is the outer bounding surface to a sheet of axisymmetric current defined by opened 
stellar field lines that reverse poloidal directions as the current sheet is crossed and 
we find ourselves in the dead zone of the overall X- wind/funnel- flow configuration 
(see Fig. 1). Until we actually construct such a boundary-layer/current-sheet theory, 
we adopt a simple modification to deal with the problem: we truncate the formal 
wind solution at some ipi < 1, below which J and (5 remain finite. We then add the 
part between ipi and 1 to the dead zone fields of the problem, i.e., treat the last few 
streamlines as opened vacuum fields and impose the pressure balance condition at 



4.3. Asymptotic Solution 

The asymptotic solution at large distances from the X-point was constructed 
by Shu et al. (1995). In particular, for a wind reaching more or less constant 
terminal velocity, its density scales roughly as p oc w~ 2 , and the Alfven discriminant 
X — > —l/f3 2 pw 2 is a slowly varying function of r. By ignoring all radial derivatives 
compared to angular derivatives, the GSE and the BE admit solutions of the form 

X = -1/PC, sin 9 = sechfCr^C, ^)], I(C, V) = jf ^ =f=^ > 

(4-10) 

where 9 is the usual polar angle in spherical coordinates, and C is a "constant of 
integration" that vary slowly in r. In an inertial frame, the wind reaches a terminal 
velocity given by 

v w = (2J -3-2Cf3) 1/2 . (4-11) 

To determine the constant C, we impose pressure balance between the X-wind and 
the dead zone. Since the dead zone field lines carry no inertia, they do not develop 
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a toroidal component, and the poloidal field satisfies the vacuum equation (Ostriker 
& Shul994). Asymptotically, we do not expect the field lines to pinch toward the 
rotational axis since the hoop stress is vanishingly small (Shu et al. 1995). For 
simplicity, we assume the boundary layer deviates only slightly from a cylindrical 
surface at the asymptotic infinity (an assumption which shall be checked a posteriori 
for consistency). For any given (large) value of r, we can approximate the boundary 
locally by w = const. Then a particular solution is B = B^ c z. For the hollow-cone 
region to trap the same amount of net flux as the wind part and have a cross-sectional 
area of 7rw^ c , we have 

R 20 h J 

-£>hc — 2 — ' 

where 0hc is a number ranging from 1 to 3 depending on the fraction of closed field 
lines in the dead zone (compare Fig. 1 of this paper with Fig. 1 of Shu et al. 2001). 
The maximal case 0hc = 3 has three times as many field lines as the minimal case 
0hc = 1) but the extra field lines cancel in oppositely directed pairs and contribute 
no net flux. The overall solution does not depend sensitively on the number 0hc, and 
we take 0hc henceforth to be unity as illustrated in Figure 1 of the current paper. 

In contrast, the wind region is dominated by the toroidal field 
R -R J ~ W2 ° 

The poloidal field B WiP = f3pv w oc 1/zu 2 is much weaker in this limit. By equating 
the magnetic pressures on both sides of the boundary, Bf lc = B 2 ^ , we obtain 

c = 2p_ = 2/3 cosh[c -i /(C) ^ (4 _ 12) 

zu hc r 

which implicitly defines C (r) . Since / only depends very weakly on C, this expression 
shows that C — > logarithmically as r — > oo. Notice that this limiting behavior 
of C ensures that wy ic deviates from a constant only logarithmically slowly, which 
validates our assumption on the geometry of the boundary layer. Written in our 
coordinates, the asymptotic geometry of each streamline is given by 

z = zusmh[C- 1 I{C,^)}. (4-13) 

In other words, each streamline is approximately radial, with a logarithmic collima- 
tion toward the axis. 

With p — > C/j3zu 2 and v w — > constant, the poloidal and toroidal Alfven speeds 
are given by 

v% P = B llP ^ C(3v 2 Jw 2 , v% v = Bl/p - Cf3. (4-14) 

Thus the Alfven speed is dominated by the toroidal component, which decreases 
logarithmically. This means that the (poloidal) terminal velocity is super-fast in 
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the asymptotic regime, and the wind has to make a fast mode transition along each 
streamline. The above analysis simply reiterates the claims made in Appendix [A] 
that the asymptotic behavior of the flow is governed by a hyperbolic differential 
equation. 



5. Global Solutions 

As a particular example, let us suppose that diffusive mass loading onto field 
lines in the X-region produces a f3 function which has the form 

P = lp(i-^r 1/3 - (5-1) 

It is easy to verify that f3dijj = f3. To be definite, let us also assume that the upper 
boundary of the X-wind near the X-point forms the maximum angle dx = 7r/3 with 
the equatorial plane in order for magnetocentrifugal acceleration to operate. The 
0[{w - 1)~ 2 ] solution (03]) takes the form 

2l = tan^ = tan-[l-(l-^) 2/3 ]. (5-2) 
3 

In fact we have chosen a very special value for the opening angle. Recall that 
the formal boundary between the X-wind and the dead zone is characterized by 
vanishing p with finite magnetic field. That results in (3 — > oo and x = ^~ 2 - If 
dx were smaller than ir/3, one may check that the boundary condition \ = w~ 2 
agrees with the series solution of §4.11 to the second order for all values of q. This 
integration constant is computed by expanding the shape of the last streamline near 
the X-point. However, when dx = vr/3, the series solution agrees with the boundary 
condition only if the quantity q in equation (14-71) satisfies 



1 (1 



7T 



If °&x > tt/3, the solution becomes discontinuous. This behavior is consistent with 
our physical intuition. When the flow is cold, the upper boundary of the X-wind is 
imposed by pressure balance. As one eases up the external pressure, dx increases. 
However, even when the external pressure drops to zero, the matter carrying stream- 
lines are confined to d < tt/3, since it is the boundary where centrifugal effects can 
overcome gravity. At least near the X-point, there is no freedom to choose the 
shape of the last streamline. Any excursion across this boundary requires additional 
pressure support from the X-wind, which calls for a warm rather than cold outflow. 

In the particular example we are studying here, the second order coefficient for 
z becomes 



z 2 = - sec 2 $ 



7 71 

12 3\/3 



(5-3) 
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For numerical tractability, we place the boundary layer at ip — 0.99. Given the 
choice of mass loading in equation (15-1 p . the (3 function is not much larger than 
unity there. 

5.1. Fixing the Free Function J(ip) 

The BE (12-13bl) is actually a quartic algebraic equation for the Alfven discrim- 
inant once the shape of the streamlines are known. The Alfven surface here is not 
a real singularity of the equation; it simply ensures that x — is a solution when 
J = w 2 . In other words smooth crossing of the Alfven surface does not uniquely 
determine the value of J. To see this, let us define 

It is always negative and asymptotes to zero for the wind since the magnetic field 
lines form a trailing spiral. The BE can be written as 

(|VVf + C 2 ){(3 2 C + J - to 2 ) 2 + 2w 2 V cS C 2 = 0. (5-4) 

As long as J is larger than some critical value, there are always real and finite 
solutions to this equation, which means the Alfven surface is automatically crossed. 
On the other hand, the fast point is a real critical point for the BE. A smooth fast 
mode transition demands the BE to have a double root at the critical point (see 
Fig|5J). That means not only does the left hand side of (15-41) need to vanish, its 
derivative with respect to £ must vanish as well. 

After some algebra, these requirements can be written as 

£ = |V^| 2/3 (J-^ 2 ) 1 / 3 r 2/3 - (5-5) 

This expression is simply a statement that at the fast point, the poloidal fluid velocity 
is equal to the magnetosonic speed, which for e = is equal to the total Alfven 
speed. Notice that both equations ()5-4p and (15-51) are automatically satisfied by 
C = (J — w 2 ) = 0. This solution, however, is unphysical since it has a discontinuity 
on the Alfven surface when x — 0. Substituting equation (I5-5P back in to the BE 
()5-4p . we have 

"|VVi 2/3 /3 4/3 + (J-^ 2 ) 2/3 

For a given value of J, the solutions to this equation give the locations where the 
BE has degenerate roots. If J < J c , equation (15-61) has no roots in the super- Alfven 
region [w 2 > J). If J > J c , then equation (I5-6P has two roots in the super- Alfven 
part of the flow. The desired solution is obtained when J = J c , and there is only 
one double root occurring at the fast mode transition point (see Fig. 2). 



2zu 2 V cS = 0. 



(5-6) 
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5.2. Interpolation Schemes and Numerical Strategy 

Our strategy is then to find interpolations between the X-point solution in 
§4.11 and the asymptotic solution of §4.31 so that the action (13-21) is extremized. 
Since the action involves an integral extending to r — > oo, and the streamlines 
are approximately radial, in general the action integral is infinite. In the X-wind 
problem, however, the assumption of stationarity is an approximation that must 
fail physically at very large distances from the X-point. If the flow extends all 
the way to spatial infinity, then steady state cannot be established in finite time. 
To make the practical aspect of this problem manageable, we opt to truncate the 
action integral at some finite spatial surface, and assume that the solution is identical 
to the asymptotic solution beyond that point. Then the interpolation requires the 
intermediate solution to join smoothly on to the asymptotic solution at the boundary. 
Since the parameter C that appears in the asymptotic solution is purely a function 
of r, it is natural to choose the boundary surface at r = tq < oo. Thus, along a given 
streamline labeled by ip, the action involves an integral over the range w G [1, ©oo], 
where 

2ff coshfC- 1 /^, 1)] 
Wco CcoshtC- 1 /^)] 1 } 

Here / is the integral defined in equation fl4-10p . and the asymptotic value of z is 
given by 

z^w , i/)) = sinh [C -1 /(C, ip)] , (5-8) 

Since the asymptotic behavior of the streamlines are predominantly radial with 
a logarithmic collimation toward the pole, we may approximate them by linear func- 
tions. There is a large class of basis functions in which z(if), w) can be expanded. 
To avoid unphysical oscillations introduced by higher order polynomial interpola- 
tions, we approximate z by a cubic spline such that the second derivative z mzu is a 
continuous piecewise linear function. 

z-ujuj — fi + (vj — w^——— — , for tUi < zu < w i+ i, (5-9) 

where i = 0...N — 1, with hj q = 1 and zun = w^. The boundary conditions on z mvo 
read 

fo = 2z 2) f N = 0. (5-10) 



Direct integration yields (j36l ) 



z = ayi + by i+1 + cfi + df i+1 , (5-11) 

Vi+i-yt 3a 2 -1 362-1 

{zu i+ i - Wi)ji H [zu. t+ i - Wi)fi + i,{b-L2) 



where 



6V — t-r± ti i j it ^ 





W i+ \ — W W — UJi 

, o = 1 — a = , 
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1 . 1 
o o 

and yi = z(wj) for that particular streamline. The y^ are determined by demanding 
z m is continuous throughout the domain. Explicitly, 

Wi - cc7 i+ i - Wi-i zu l+ i —Wi y i+ i - yi yi - . . 
7 + o fi+ a /<+! = ' ( 5 " 13 ) 

which is a set of iV — 2 linear equations for the N y^ The boundary conditions yo — 
and T/jv — -Zoo close the equations, and allows unique determination of y^ once /j are 
given. Since we have information on the slope of the solution on both boundaries, 
they impose two further constraints 

Zl = M — _ h Wl _ Wo )/ -h Wl - O7 )/i, (5-14) 

V3\ — Wq 6 O 

Vn — VN-1 1/ w . 1/ w /_ 1K s 

^oo,w = 1- tI^jv - ■cjn-i)jn-i + 77 (^jv - ^N-IJJN- (5-15) 

To demonstrate the principles, we choose N = 3, so that all the /j are constrained. 
For a given set of the equations (I5-I3p . (15-141) . and (15-151) form a set of four 
linear equations, which can be solved by standard means. We also define ©2 by 

T = ' ( 5_16 ) 

TJJ\ — 1 "072 ~~ ^1 

i.e., we demand that the interval between interpolation points to increase expo- 
nentially. Thus the shape of each streamline is parameterized by a single variable, 

The action integral and the asymptotic solution can be treated as solutions to 
a set of simultaneous "ordinary" differential equations 

ds r°° w di p 

— = Lwaw, — = — . =, (5-17) 

dip Ji dtp y/2J(7p)-3-2CP(ifj) 

subject to the boundary conditions 

S(0) = 0, 1(0) = 0. 

Here L represents the Lagrangian appearing in the action (13-2 j) . For each value of 
ip, to compute the right hand side of equation (I5-I7p . one needs the values of Wi(ip), 
I(ip) and I(ipi), where ip\ = 0.99 is the label of the boundary layer discussed in 
§4.21 Ideally, one would like to specify the shape of the last streamline by fixing the 
values of wi(ipi) and I(tpi) as boundary conditions, and vary the function vj\(ijj) in 
a constrained manner to achieve a local extremum of the action. In practice, we find 
it more convenient to implement a scheme where only c^i^i) is given, and I(ipi) is 
determined as an eigenvalue. This approach allows more freedom in the parameter 
search for the desired wi(ip). With each streamline fully parameterized, one can 
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proceed to determine the necessary value of J(if>) that allows a smooth fast mode 
transition according to the procedure outlined in §5.11 



Once J (if}) and I(ip) are both known, we can easily solve the BE (12-13bj) as 
an algebraic equation along each streamline for £ using standard techniques such 
as Laguerre's method (see Press et al. 1992). In particular, note that we do not 
actually use the extremal property of the action principle with respect to A to 
attack the Bernoulli equation, but effect direct solutions of it instead. Increased 
numerical accuracy constitutes only one reason for a mixed procedure, where we do 
find the extremal action through variations of if), or equivalently, through variations 
of z(zcr,if)), as a substitute for solving the Grad-Shafranov equation. There is a yet 
more practical reason. It turns out the the correct solution sits on a saddle, where 
the extremal action is minimized by variations of if> but maximized by variations of 
A. This combination makes a numerical search for the extremal action extremely 
difficult to execute in practice, perhaps even impossible, if the search is carried out 
in the double-function space of allowable if) and A. 

One further obstacle to overcome is that the action integral (!3-2p is logarithmi- 
cally divergent at the X-point. Recall that the Jacobian of the coordinate transfor- 
mation vanishes at the X-point since it maps the entire axis of w = 1 onto a single 
point. A series expansion of the Lagrangian using the series solution of section §4.11 
shows that it diverges as 

L - = v£*ry < M8 > 

Fortunately, this term does not enter into the variation scheme, and we may safely 
remove it as a counter term from the Lagrangian, as is the standard practice in 
quantum-field theory. 



Finally, the function Wi(ip) is modeled by a Hyman filtered spline (1211 ) inter- 
polating over evenly spaced control points if>i G [O,^]. The values of W\ at these 
control points, Wi(if>i), are the parameters we can adjust in our variation scheme. We 
restrict the parameter space to that satisfies the condition that the streamlines do 
not cross and that each streamline is monotonic. We then adopt a genetic algorithm 
to search for a set of wi(if>i) that gives a local extremum of the action (13-21) . 



6. Numerical Results 

We compute the streamlines for three cases of average mass loading correspond- 
ing to /3 = 1, 2, 3. In each case, we place the outer boundary of the computational 
domain at a constant radius so that it intersects the last streamline at vj^ipi) = 20 
(which yields C = 0.1(3). After a multidimensional search, we locate the desired 
set of control points that extremize the action. They are tabulated in Table HJ an d 
the function zui(if>) is interpolated between these points as described in the previous 
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section. 



p 


O7i(0.0) 


O7i(0.2) 


roi(0.4) 


O7i(0.6) 


O7i(0.8) 


n7i(0.99) 


1.0 


29.687 


29.817 


23.281 


18.809 


8.310 


6.000 


2.0 


28.281 


29.165 


24.644 


17.774 


11.150 


6.000 


3.0 


28.384 


28.985 


23.990 


19.965 


10.139 


6.000 



Table 1: Values of control points W\{ipi) that yield a local extremum of the action. 
The last value Wi(0.99) is fixed as a boundary condition. 



For each converged solution, we can numerically integrate the asymptotic equa- 
tion to evaluate I(ip). For practical purposes, we present here an interpolation 
formula that is a seventh degree polynomial in and the coefficients are tabu- 
lated in Table [2J Once is known, one may determine the outer boundary of the 
computational domain in accordance with the asymptotic condition (14-101) . With 
the combination of wi(ip) and I (if)), we are able to reconstruct the streamlines with 
the spline interpolation scheme, and they are depicted in Figure [3j The location of 
the Alfven surface determines the value of J as a function of ip, which ultimately 
allows us to compute the angular momentum being transported as well as the termi- 
nal velocity along each streamline. For convenience, we also present an interpolation 
formula for J(ip) as a polynomial in /3, with the coefficients tabulated in Table [3j 






io 


h 


h 


h 


h 


h 


h 


h 


1 


0.732 


-0.446 


0.886 


-0.511 


-1.781 


2.648 


-1.397 


0.263 


2 


0.842 


0.413 


-2.504 


2.699 


-11.768 


24.849 


-22.870 


7.602 


3 


1.164 


-2.915 


30.674 


-194.996 


585.524 


-998.495 


930.236 


-347.764 



Table 2: Interpolation formula for Ii nt (ip) = ^2 i=Q Iif3 % - The interpolated function 
agrees with the numerical values to within 0.5%. 






Jo 


h 


J 2 




Ja 


h 




Jr 


J 


1 


-2.791 


17.214 


-13.640 


-10.294 


22.483 


-13.656 


3.628 


-0.362 


2.638 


2 


-14.944 


45.162 


-43.065 


21.660 


-6.265 


1.057 


-0.0969 


0.00373 


4.356 


3 


-20.285 


40.676 


-25.088 


8.006 


-1.440 


0.148 


-0.00825 


0.000191 


6.202 



Table 3: Interpolation formula for Ji n t(^) = Y^i=oJi0 l - The last column gives the 
value J of J{ip) averaged over ip from to 0.99. The interpolated function agrees 
with the numerical values to within 1%. 

The solid lines in Figure 3 show the logarithmically spaced contours of constant 
density. It is evident that even though the dotted streamlines become asymptotically 
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radial and only collimate logarithmically slowly, the density becomes cylindrically 
stratified very quickly, giving the X-wind the illusion of a jet-like appearance (Shang 
et al. 1998, 2002). 

Detailed comparisons of the results obtained here with those given by Shang 
(1998) show some differences, but the main impression is how remarkably well the so- 
lutions obtained by the two very different methods for the same mass-to-flux loading 
f3{ip) agree with one another. Shang (1998) had a similar experience in comparing 
her approximate, but analytic, solutions for the sub-Alfvenic region to the exact, 
but numerical, solutions obtained by Najita & Shu (1994). 

We attribute the fortunate circumstance to the following causes. If one is given 
somehow the geometric shape of the streamlines (or, equivalently, the field lines 
in the meridional plane), then the Weber-Davis procedure used by Shang, which 
includes an exact solution of Bernoulli's equation, would give an exact solution of 
the two-dimensional flow problem, provided one takes care to cross each of the critical 
points properly. In realistic circumstances, the geometric shape of streamlines in the 
meridional plane is not given a priori, but is to be found from the Grad-Shafranov 
equation (or, equivalently, from minimizing the action by variations of the stream 
function ip). However, if one has analytic solutions to the Grad-Shafranov equation 
(from the work of Shu et al. 1994b and 1995) near and far from the X-point, then 
there are only so many ways that one can adjust the function z(zu,ip) for values of 
ip from to 1 and of w close to 1 (or dimensionally, R x ) to > 1 (or R x ) that 
will connect the shape of the streamlines near the X-point (a fan) smoothly to those 
appropriate at asymptotic infinity (radial outflow). The procedures used by Shang 
(1998) and those used here to make such adjustments differ, but the global solution 
is relatively insensitive to these details as long as one gets the conserved quantities: 
mass-to-flux loading (3(i/j), angular momentum distribution J(ip), and Bernoulli's 
constant H(ip) = correctly. 



For the convenience of the reader, we summarize the recipes needed to convert 
the results of the previous section into numerical X-wind models for astronomical and 
meteoritical applications. Begin with the equation that describes the dimensionless 
locus of a streamline for given ip with numerical value between and 1: 



7. 



Discussion and Conclusions 



7.1. Recipe for Use of Results 



Z = z(w, '0), 



(7-1) 



where the functional form of z(zu, ip) is computed numerically by the technique 
described in §5.2. 
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The reconstruction of streamline shapes, i.e., the function z(zcr, ip), is performed 
over three radial intervals whose end points are vjq = 1, voi(ip) > 1, cc7 2 (V ; ) > wi^), 
and &73(ip) = w^tjj) > W2{ip) that give a geometrically increasing separation: 

CC7i — 1 — ro l 

where zu^ip) is given by equation (15-71) : 

2^ coshjC- 1 /^,!)] 

Ccosh^- 1 /^,^)]' 1 J 

For practical computations, we choose C = 0.1/3 so that = 20 on the ■0=1 
streamline. The asymptotic integral J(C, ^) in equation (14-101) can be approximated 
by a seventh degree polynomial in (3 . 

/(V) = Jo + Ji/T 1 W + • • • + Jr/T 7 W, (7-4) 

where the coefficients Jo, Ji, . . ., 1-j are given in Table 2 for the three values of 
(3 = 1, 2, 3. The function wi(ip) represents the first nontrivial abscissa of the 
spline beyond the X-point for each value of ip and is tabulated in Table [1] for ipi = 
0.0, 0.2, 0.4, 0.6, 0.8, 0.99. For intermediate values, we interpolate zui by a piecewise 
cubic polynomial: 

m = h (ipi) + hi(ip)(i) - if)i) + h 2 {ipi)(il) - ipij 2 + h 3 (4>i)(ij - ^) 3 fon/>i < ip < ^i+i- 

(7-5) 

In Table HI we list the values of hj(ipi) for each case of (3. To get vj2(ip) for any value 
of ip, one should use equation (17-21) after first computing Wi(ip) and ■uj 00 {ip) at the 
desired value of ip. 

The shape of each streamline given by ip = const in the three radial intervals 
whose end points are wo(ip) = 1, zui(tp), zn^ip), and U7 3 (ip) = Woo(^) is then 
described by a piecewise cubic polynomial, whose form, suppressing the implicit 
dependence on ip, is given by equation (15- lip : 

z{w) = Vl a + y 2 b + fo+i-* 7 ') 2 [ /l(a 3 _ a) + /2 ( 6 3 _ 6) ] ; (7 . 6) 

where 

a = , b = . (7-7) 

The coefficients y±, y 2 , fi, and f 2 are listed in Table [5]for discrete values of ip = 0.0, 
0.2. , 0.4, 0.6, 0.8, 0.99 in the three cases j3 = 1, 2, 3. A Hyman limited spline may 
be used to compute the streamlines for other values of ip. 

The partial derivatives of ip with w or z are now given by the usual rules of 
multivariate calculus: 

df\ __{dz/dw) ± fdf\ i 

dw) z {dz/d^)J \dz) w {dz/d^) m { ' } 
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0.0 


0.2 


0.4 


0.6 


0.8 


p = 


1 




29.687 


29.817 


23.281 


18.809 


8.310 








17.315 


-1.950 


-27.520 


-37.428 


-31.809 








-153.650 


-333.100 


126.938 


-254.103 


103.428 








351.625 


897.250 


-505.688 


893.830 


0.000 


(3 = 


2 


/i 


28.281 


29.165 


24.644 


17.774 


11.150 








17.933 


-9.093 


-28.478 


-33.735 


-30.035 






h 


-67.563 


-105.763 


-61.800 


-9.272 


15.422 






h 3 


0.000 


191.000 


162.188 


61.737 


0.000 


P = 


1 


h 


28.384 


28.985 


23.990 


19.965 


10.139 






hi 


16.995 


-9.015 


-22.550 


-34.628 


-35.107 






h 2 


-79.800 


-171.725 


96.763 


-215.142 


70.117 






h 


49.250 


459.625 


-423.188 


713.150 


0.000 



Table 4: Interpolation coefficients for W\{ip) 



Table 3 gives J(ip) as a seventh order polynomial in /3(ip): 

J(V0 = Jo + JM) + ■■■ + J 7 P 7 W, (7-9) 

where (3{ip) is itself given by 

/?(V0 = !j9(1-^)- 1/s , (7-10) 

with [3 — 1, 2, 3 in the three chosen model cases. The coefficients tabulated in Table 
3 give a J(ip) that guarantees that equation ( 12-I3bj) . 



V-Ju)' =0 - (7 ' u) 

has one real root for A in the computational domain when V e s is given by equation 
(E2D: 

By solving equation (17- lip as a fourth-order polynomial, we may obtain the relevant 
value for the Alfven discriminant A. Then the density can be computed through 
equation (12-9 j) 

p = (f3 2 - w 2 A)-\ (7-13) 
Note that this equation produces p = (3~ 2 at the Alfvenic transition ^4 = 0. 

With the density in place, we may obtain the two components of dimensionless 
poloidal velocity from the definition ( 12-31) of ip: 

_ I dip _ 1 dip 
w ro = — n-, u z = — . (7-14) 

■cup oz wp aw 
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T 


0.0 


0.2 


0.4 


0.6 


0.8 


0.99 


= 


1 


it i 


0.0 


12.764 


21.872 


41.255 


51.367 


525.487 








0.0 


366.579 


627.01 


1044.15 


1402.6 


3648.96 








0.0 


3.834 x 10" 3 


3.083 x 10~ 2 


0.256 


4.120 


118.485 






fa 


0.0 


-2.901 x 10~ 5 


-4.685 x 10" 4 


-1.033 x 10~ 2 


-0.301 


-48.051 


P = 


2 


y i 


0.0 


12.091 


22.694 


29.659 


37.935 


73.230 






2/2 


0.0 


86.922 


156.801 


208.943 


267.766 


383.996 






/i 


0.0 


2.359 x lO- 3 


1.489 x 10~ 2 


9.705 x 10~ 2 


0.750 


9.943 






/a 


0.0 


-1.991 x 10" 4 


-1.518 x 10~ 3 


-1.333 x 10~ 2 


-0.158 


-4.246 


P = 


3 




0.0 


14.138 


22.415 


33.841 


26.286 


44.094 






2/2 


0.0 


82.541 


104.01 


133.299 


131.469 


173.659 






/i 


0.0 


1.903 x 10~ 2 


2.353 x 10~ 2 


5.056 x 10~ 2 


0.376 


2.950 






/a 


0.0 


-3.651 x 10~ 3 


-5.610 x 10~ 3 


-1.722 x 10~ 2 


-0.102 


-1.424 



Table 5: Spline coefficients for the streamlines. 



The toroidal velocity in the corotating frame is given by equation (12-6 j) 

Note that J(ip) = vj 2 where (3 2 p = 1 keeps the toroidal velocity u v well-behaved 
across the Alfven surface, which is not one of the critical surfaces of the overall 
problem. 

The vector magnetic field may now be obtained from equation f!2-4p : 

B = Ppu, (7-16) 

whereas the azimuthal velocity in the inertial frame is given by 

v !fi = u !fi + w, (7-17) 

with the term w from the frame rotation being cancelled at large w where — > —zu 
because p vanishes as 1/w 2 at large distances from the rotation axis. Finally, 
to convert the computed quantities to their dimensional counterparts, we must 
multiply velocities, densities, and magnetic fields by Rx^x, M w /4rrR x Qx, and 
(Q x Mw/Rx) 1/2 , respectively. 

For interpolations or extrapolations in /3, we recommend computation first of 
the dimensionless density, velocity, and magnetic fields for the three cases (3 = 1, 2, 
3, and then direct interpolations or extrapolations of those fields. Other techniques 
starting farther back in the process run the danger of obtaining complex roots of A 
(i.e., complex values of p) from the solution of the quartic equation (17-111) because 
of slight inaccuracies in computing the numerical coefficients. 
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7.2. Summary 



In this paper, we have presented a technique by which solutions to the so-called 
Grad-Shafranov equation for X-wind flow can be solved, not by attacking the par- 
tial differential equation directly, but by choosing trial functions that minimize an 
appropriate action integral. While this method has been applied before in problems 
of plasma confinement in the fusion community, we believe that the example given 
here is its first application in astrophysics for the notorious case when magnetohy- 
drodynamical flows cross critical surfaces that change the character of the underlying 



Many empirical arguments suggest that funnel flows and X-winds do underlie 
the accretion hot-spots, jets, and winds of YSOs, although a dipolar field geometry 
near the star (see Fig. 1) may be an over-simplification (Ardila et al. 2002, Unruh et 
al. 2004, Johns-Krull 2007). Fortunately, although the fractional areal coverage of 
hot spots depends on the detailed multipole structure of the surfaces of actual young 
stars, the general validity of X-wind theory depends only on the level of trapped flux 
in the X-region and is insensitive to the magnetic geometry on the star as long as 
the fields are strong (Mohanty & Shu 2007). The trapped flux in the X-wind models 
of this paper are computed as 



and should be compared with the magnetic flux (area times mean field) in hot-spots 
on one hemisphere's surface of the star impacted by the corresponding funnel flow. 
(Both fluxes are 1/3 of the total trapped flux in the X-region and equal the net flux 
of the dead region.) For T Tauri stars, the comparison is pretty good (see, e.g., 
Johns-Krull & Gafford 2002). 

Apart from relative simplicity, the semi-analytical solutions summarized in §7.1 
have many other advantages. For example, the solutions hold over a formally infi- 
nite dynamic range, showing the asymptotic, logarithmically slow, collimation into 
jets missing in many numerical simulations. These properties make the models of 
this paper especially suitable for a wide variety of astronomical and meteoritical 
applications, such as detailed comparisons with observations, trajectories of solids 
entrained in the wind, and interactions with neighboring circumstellar or interstellar 
matter. A needed generalization for future research is the inclusion of the effects of 
the intrinsic magnetization of the surrounding accretion disk. 



PDE. 




(7-18) 
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A. Character of Governing Equation 



The GSE (12- 7p resembles the steady state heat diffusion equation with a variable 
diffusion coefficient A. This analogy is actually misleading since we do not know its 
overall character until we substitute in the implicit dependence of A on ip by solving 
the (algebraic) BE and examine the characteristics of the GSE. To do so, let us first 
differentiate the BE with respect to w and z. 



.zvd ) 

2roVU 



2A, 



J 



((3 2 - w 2 Af 

,zz, 

2w A A, z 
+ {(5 2 -w 2 Af 



A* 

2V eS + 2e 2 In 
2A 



e 2 h 



A 3 



(3 2 - w 2 A 

2 



+ 



2V eS + 2e 2 In 



e 2 h 



p 2 - w 2 A 



+ 



0. 



where in the above equations, a subscript denotes partial derivative and the ellipsis 
symbols include terms that are irrelevant in determining the character of the GSE. 
These equations may be solved for A, m and A z to give 



A 



A, 



where 



V 



zu 



[3 2 - w 2 A 



A 3 V^ 2 J 



P 2 



e 2 w 4 



(3 2 - w 2 A ((3 2 - w 2 Af 



after we eliminate V e s in the expression by using the BE (12-I3bl) . The second 
derivative terms in the GSE (I2-I3al) can now be written in the form 



a^,mvo + 2bip <u?z + c^, zz + ... = 0, 



where 



A + 



C = A + 



j> 1 V V 

The character of the GSE is determined by the quantity A = b 2 — ac (Garabedian 
1986): it is elliptic, parabolic, or hyperbolic if A is negative, zero, or positive. We 
may compute A for our GSE explicitly. 



A 



-A 2 



+ (Jw- 



1) 2 A- 



e 2 Aw A [(3{(3 2 -w 2 A)\ 



w 2 A(3~ 2 \Vi)\ 2 + (Jw- 2 - 1) 2 A~ 2 - e 2 Aw i [(3((3 2 - w 2 A)\ 



■ (Al) 
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The interpretation of this expression becomes transparent if we transform back into 
the physical quantities. After some algebra, we have 



A = 

where va p 
with B 2 = 



-A 2 



U 2 -6 2 (l-M%) 
L(l-Mi)(^- e 2) +M 2j 



A 2 



_(uj-vl p ){uj-vlp)_ 



(A2) 



B 2 /p is the poloidal component of the Alfven velocity va = a/ B 2 / p 
Up is the poloidal fluid velocity, and v s is defined by 



2„,2 



e v 



Ap 



In the limit where v 2 A 3> e 2 , it reduces to the thermal sound speed. In addition, v± p 
denote the poloidal component of the fast and slow MHD wave speeds, respective, 
and are given by 



l(v 2 A +e 2 



1 ± 



4w 2 



(A3) 



A moment of thought reveals that v, 
flA2p is now clear. The governing GSE is elliptic when u 2 
and it is hyperbolic when v 2 < u 2 < v\ 



< V-p < v +p . The significance of the equation 

^< v 2 or v 2 _ p <u 2 p < v 2 +p , 



[ p or u 2 p >v 2 +p dl9|;l39|). 



To be definite, we shall refer to the loci where the poloidal velocity squared 
equals v 2 , f 2 ^ and v\ p as sonic, slow, and fast surfaces, respectively. Despite the 
deceiving appearance of the GSE (12-71) . note that it does not change character on the 
Alfven surface when ^4 = (or equivalently when Ma = 1 and the total fluid velocity 
in the corotating frame equals the total Alfven speed); it remains elliptic until the 
fast surface. The fact that the asymptotic flow is described by a hyperbolic PDE is 
consistent with our physical intuition. When the fluid speed is super-magnetosonic, 
no information can be sent upstream into the flow. Thus the asymptotic behavior 
of the X-wind is determined by the "initial condition" at the place when the fluid 
velocity first becomes equal to the fastest signal propagation speed - a defining 
feature of hyperbolic problems. 



In the cold limit the discriminant A has the simplification, 



A 



-A 2 



|VVf + (Jw 



-2 



l?A- 



w 2 A(3~ 2 \Vi)\ 2 + (Jw- 2 - IfA- 



-A 2 



u 2 p /v 2 A 



(A4) 



This equation explicitly states that the transition to the hyperbolic portion of the 
solution is done through the fast surface, where the poloidal fluid velocity is equal to 
the magnetosonic speed, which is the total Alfven speed va = B / \fp when e is set to 
zero. The axial symmetry of the assumed problem guarantees that any compressions 
or rarefactions occur only in the meridional plane, so the relevant speed of signal 
propagation in the limit e — > is the magnetosonic speed relative to the poloidal 
motion of the fluid. 
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Fig. 1. — Funnel flow (red curves), X-wind streamlines (blue curves), and field lines 
dead to magnetocentrifugal fling (black curves) according to Ostriker & Shu (1995) 
and Shang et al. (1998). The magnetic field near the origin (center of YSO) is 
modeled as a magnetic dipole, and all of the field lines contained in the X-wind 
have their counterparts (with reversed directions) in opened stellar field lines that 
lie inside a hollow cone dead to flow surrounding the z-axis. Exact pressure balance 
across the sheet current that divides the X-wind and dead field lines holds near and 
far from the Y-point at w ~ 1.3, z ~ 0.7, but this balance is only approximate 
at intermediate distances, which accounts for why the tilted upside-down Y of the 
separatrix does not have the equal angles of 120° that would characterize an exact 
Y-configuration appropriate for coronal conditions. In fact, the field lines in red 
and black are computed as if they were vacuum magnetic fields; only the portion 
depected in blue are attacked via the solution of the Grad-Shafranov equation or 
its variational-principle analog discussed in the present paper. The neutral line 
separating the black and blue field lines is replaced by a separatrix of prescribed 
locus satisfying the approximate pressure balance as described above, (see §4). 
The implied poloidal and toroidal current flows are discussed in Ostriker & Shu 
(1995) and Shu et al. (1995); in particular, there is no current flow along the z-axis 
which is taken to be a region of vacuum longitudinal field. 




Fig. 2.- 



Determination of the critical value of J that allows a fast mode transition 
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Fig. 3. — Solutions for (3 = 1,2,3. The dotted curves represent the streamlines 
labeled by constant ip, and the solid curves are isodensity contours separated by 



